#!/usr/bin/python3
# -*- coding: utf-8 -*-
"""
@FileName: Upper_Hessenberg.py
@Version: 1.0
@Date: 2023/7/14 10:20
@Description:
"""

import numpy as np
import numpy.linalg as npla
import math


def hessen(a_: np.ndarray):
    eps = 1E-15
    m, n = a_.shape
    # 每一列是 HouseHolder 矩阵的生成向量
    v_ = np.zeros((m, m))
    # 主对角线 m 个元素; 副对角线 m-1 个元素, 副对角从第二行开始, python 索引为1
    # 最后移动的列向量是 2x1 维; 所以截止行索引是 m-2, python 索引是 m-1 (<右索引)
    for k in range(1, m - 1):
        # 取副对角线行(-1) 往下的整列
        x = a_[k:m, k - 1:k]
        # 第一行长度为 m-1, 第二行为 m-2, ...; 只有首元素不为零; sign 是为了避免临近数字相减
        v_tmp = -np.sign(x[0] + eps) * npla.norm(x) * np.eye(m - k, 1) - x
        #  归一化, 则 HouseHolder=I-2P, P的分母=1
        v_tmp = v_tmp / npla.norm(v_tmp)
        # 左乘 H 矩阵: 上面 [0:k,0:m] 的子块不变; 下面 [k:m,1:k-1] 的子块也不变, 左侧元素已被移动为零
        a_[k:m, k - 1:m] = a_[k:m, k - 1:m] - 2 * v_tmp @ v_tmp.T @ a_[k:m, k - 1:m]
        # 右乘 H 矩阵: 左边[0:m,0:k]的子块不变
        a_[0:m, k:m] = a_[0:m, k:m] - 2 * a_[:, k:m] @ v_tmp @ v_tmp.T
        # record H矩阵的生成向量
        v_[0:m - k, k - 1:k] = v_tmp
    return a_, v_


if __name__ == '__main__':
    a = np.array(
        [
            [4, 3, 1, 1],
            [2, 6, 2, 1],
            [1, 4, 5, 8],
            [8, 5, 6, 9]
        ]
    )
    a, v = hessen(a)
    print(f'Upper hessen:\n{a}')
    print(f'house holder:\n{v}')
